#Salvar a base de dados (BD_survive - recortada_final.csv) em uma pasta no coputador
#Ao usar o comando, read.csv(file.choose(), header = TRUE,sep=";", dec = ","),
#abrirá uma janela para o pesquisador indicar em que pasta do seu computador
#está salva a base de dados
dados<-read.csv(file.choose(), header = TRUE,sep=";", dec = ",")

dim(dados)
head(dados)

names(dados)
tempo<-dados$tdias
censura<-dados$event
require(survival)

#Estimacao da Funcao de sobrevivencia pelo método de Kaplan-Meier
KM<-survfit(Surv(tempo,censura)~1,conf.int=F)
summary(KM)
plot(KM,conf.int=F, xlab="Tempo", ylab="S(t)",mark.time=T)

#Estudo do comportamento da funcao de risco dos tempos
require(AdequacyModel)
TTT(tempo)

#Estudo de qual(is) modelo(s) parametrico(s) se ajusta(m) bem aos dados
##### 1º) Distribuicao exponencail  ###
#Estimação dos parâmetro da Distribuicao
mexp<-survreg(Surv(tempo,censura)~1, dist='exponential')
mexp

alphaexp<-exp(mexp$coefficients[1])
alphaexp

n<-length(tempo)
pexp<-1

#Calculo das medidas AIC, AICcorrigido e BIC
AICexp<-(-2*mexp$loglik[1])+(2*pexp)  

AICcexp<-AICexp + ((2*pexp*(pexp+1))/(n-pexp-1))

BICexp<-(-2*mexp$loglik[1]) + pexp*log(n)

medidaexp<-cbind(AICexp,AICcexp,BICexp)
medidaexp


##### 2º) Distribuicao Weibull ###
#Estimação dos parâmetro da Distribuicao
mwe<-survreg(Surv(tempo,censura)~1, dist='weibull')
mwe

alphaw<-exp(mwe$coefficients[1])
alphaw

gama<- 1/mwe$scale
gama
pws<-2

#Calculo das medidas AIC, AICcorrigido e BIC
AICws<-(-2*mwe$loglik[1])+(2*pws)

AICcws<-AICws + ((2*pws*(pws+1))/(n-pws-1))

BICws<-(-2*mwe$loglik[1]) + pws*log(n)

medidasw<-cbind(AICws,AICcws,BICws)
medidasw

## Teste da Razao de Verossimilhanca###
#Comparando as distribuicoes Exponencial e Weibull
TRV<-2*(mwe$loglik[1]-mexp$loglik[1])
1-pchisq(TRV,1)

##### Estimativas da sobrevivencia usando Kaplan-Meier, modelo exponencial, Weibull#####

km<-survfit(Surv(tempo,censura)~1)
time<-km$time         ### tempo de KM ###
#time

skm<-km$surv     ### funcao de sobrevivencia estimada de Kaplan-Meier ###
#skm

sexp<-exp(-(time/alphaexp))  ### funcao de sobrevivencia estimada da distribuicao exponencial ###
#sexp

swe<-exp(-(time/alphaw)^gama)  ### funcao de sobrevivencia da distribuicao weibull ###
#swe

#Graficos das funcoes de sobrevivencia estimadas acima para comparar
plot(km,conf.int=F, xlab="Tempos", ylab="S(t)",mark.time=T)
lines(c(0,time),c(1,sexp),lty=2,col=2)
legend(800,0.9,lty=c(1,2),col=c(1,2),c("Kaplan-Meier","Exponencial"),bty="n",cex=0.8)
lines(c(0,time),c(1,swe),lty=3,col=3)
legend(800,0.8,lty=c(3),col=3,c("Weibull"),bty="n",cex=0.8)

#Ao considerar as analises, gráfica e medidas (AIC, AICc, BIC), a conclusao é que o modelo Weibull
#se ajusta melhor aos dados

#grafico apenas com ajuste da distribuicao Weibull
plot(km,conf.int=F, xlab="Tempos", ylab="S(t)",mark.time=T)
lines(c(0,time),c(1,swe),lty=3,col=3)
legend(800,0.9,lty=c(1,3),col=c(1,3),c("Kaplan-Meier","Weibull"),bty="n",cex=0.8)

#Residuo de Cox-Snell para verificar a qualidade global do ajuste do modelo weibull#
swe<-exp(-(tempo/alphaw)^gama)
ei<-(-log(swe)) # residuo de Cox-Snell

KMew<-survfit(Surv(ei,censura)~1,conf.int=F)
te<-KMew$time  # residuo de Cox-Snell
ste<-KMew$surv
sexp<-exp(-te)

plot(ste,sexp, xlab="S(ei):Kaplan-Meier", ylab="S(ei):Exponencial padrao")
plot(KMew,conf.int=F, xlab="Residuos de Cox-Snell", ylab="Sobrevivencia estimada",mark.time=T)
lines(te,sexp,lty=2, col=2)
legend(0.6,0.8,lty=c(1,2), col=c(1,2),c("Kaplan-Meier", "Exponencial padrao"), cex=0.8, bty="n")

#O Residuo de Cox-Snell confirma a adequabilidade do uso do modelo Weibull
#para analisar os dados em estudo


